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Abstract. In numerical studies of the dynamics of unbound quantum mechanical 
systems, absorbing boundary conditions are frequently applied. Although this 
certainly provides a useful tool in facilitating the description of the system, its 
applications to systems consisting of more than one particle is problematic. This 
is due to the fact that all information about the system is lost upon absorption of 
one particle; a formalism based solely on the Scrhodinger equation is not able to 
describe the remainder of the system as particles are lost. Here we demonstrate how 
the dynamics of a quantum system with a given number of identical fermions may 
be described in a manner which allows for particle loss. A consistent formalism 
which incorporates the evolution of sub-systems with a reduced number of particles 
is constructed through the Lindblad equation. Specifically, the transition from an N- 
particle system to an (N — 1) -particle system due to a complex absorbing potential 
is achieved by relating the Lindblad operators to annihilation operators. The method 
allows for a straight forward interpretation of how many constituent particles have 
left the system after interaction. We illustrate the formalism using one-dimensional 
two-particle model problems. 



PACS numbers: 31.15.-p,03.65.Ca,32.80.-t 
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1. Introduction and basic theory 

Large efforts have been invested in understanding the quantum mechanical behaviour 
of dynamical unbound systems involving several particles. Experimental advances 
allow us to study processes in which internal interactions between the constituent 
particles play a crucial role - in addition to dynamics induced by external 
perturbations. Examples of such processes may be photo ionization of helium HI and 
cascaded Auger processes following ionization of inner core electrons in atoms 0. 
Of course, in order to understand such processes, theoretical and numerical studies are 
required. Furthermore, the theoretical interest of such systems, be it within the context 
of solid state, molecular, atomic or nuclear physics, is spurred by the fact that they 
represent demanding tasks. One obvious challenge in describing unbound systems is 
that their spatial extension may become arbitrarily large. Moreover, even if one is able 
to represent the whole system within a finite space, extracting relevant information 
from the final wave function may be far from trivial. 

The unbound quantum systems under study may often be thought of as having 
an interaction region of finite spatial extension and an asymptotic region where the 
unbound part of the system travels outwards. As this asymptotic behaviour often 
is well known, it may be desirable to describe only the dynamics of the part of the 
wave function belonging to the interaction region. In a numerical implementation this 
cannot be done by simply resorting to a representation of space that is smaller than the 
extension of the wave function as this leads to unphysical reflections at the boundary. 
However, it can be achieved by imposing absorbing boundary conditions, i.e. by 
demanding that the wave function vanishes as the particle approaches the edge of the 
numerical grid - preferably with as little reflection as possible, see references (31 31 
or the recent review of reference 0. Such absorbers are frequently referred to as 
perfectly matched layers in the context of general wave equations (6l. 

When propagating a wave packet on a numerical grid, a common way to impose 
absorbing boundary conditions is by adding a complex absorbing potential (CAP) to 
the Hamiltonian of the system [0. Complex absorbing potentials are widely used e.g. 
within molecular dynamics (Hill CGI ant ^ atom i c physics [TTTl[T2l . Alternatively, this 
way of absorbing particles may be formulated as multiplying the wave function with a 
masking function at each time step fl}. 

In any case the resulting effective Hamiltonian acquires an anti-Hermitian part, 



where both H and T are Hermitian and T is positive semi-definite (T > 0). 
The wave function of the system obeys the Schrodinger equation 



H eff = H-iT, 



(1) 



ih%-Mt))=H eS \V{t)). 



(2) 




(3) 
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It is easy to see that the evolution is non-unitary, viz, 

i T r^(0] = -^Tr[r P (0]<0, (4) 

so that probability is "lost" if the density matrix overlaps with anti-Hermitian part of 
the effective Hamiltonian. 

Although methods involving non-Hermitian effective Hamiltonians may come in 
very handy in reducing the complexity in describing potentially unbound systems, 
there is one major problem when the system contains more than one particle: As one 
particle leaves the rest of the system and is subsequently absorbed, the entire wave 
function is lost. This is obvious since the wave function is normalized to the probability 
of finding all particles within the space defined by the numerical implementation. In 
other words: If an initial A^-particle system gradually "loses" one particle through some 
non-Hermitian "interaction", the corresponding wave function gradually goes 

to zero - not to some wave-function corresponding to N — 1 particles. 

As it is desirable to be able to describe dynamics where several particles are lost 
one by one, one may try and construct a formalism where an (N — 1) -particle wave 
function (t)) is created as one particle is lost, and where this wave function may 

be propagated using the corresponding Schrodinger equation including some source 
term. Once being able to do this, the extension to N — 2 particles etc. follows by 
induction. However, as it turns out, such a construction is not possible due to the fact 
that the process of losing particles in this way is irreversible: As a particle is absorbed, 
information is irretrievably lost. Hence, a Markovian master equation should be a more 
suitable starting point than a pure state approach. Lindblad 0"3l was able to show that 
in order for the evolution of a system to be Markovian, trace-conserving and positive, 
the density operator p(t) has to obey an equation of the form 

ihj f p = [H,p] - ®{p) (5) 
with the so-called Lindbladian *3>{p) given by the generic expression lfl4l 

^(P) = i Y, V m > n ( A li A nP + P A m A n ~ ^ A nP A Vj ■ (6) 

Here, H is the Hamiltonian governing the unitary part of the evolution. The operators 
A n are referred to as Lindblad operators. In a diagonal representation, the Lindbladian 
simplifies to 

®{p) = UlA n p + P AlA n - 2A nP Aj) . (7) 

n 

Equation © is usually used to describe energy dissipation to the environment 03). 
However, it has recently been demonstrated that the spontaneous decay of unstable 
particles may also be described within this formalism |[T6l [T71 . In reference lTT6) a 
master equation of Lindblad from is obtained for a Hilbert space consisting of unstable 
particle states and vacuum, however without any dynamical degrees of freedom. The 
decay rates of the unstable particles are introduced as parameters. Similarly, in 
reference [fT7ll decay from one system to another one consisting of its decay products is 
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described by introducing a single Lindblad operator. In these works it is demonstrated 
that one may describe decoherence along with decay in this context. 

In the present work these ideas have been used to generalize the standard one- 
particle absorbing boundary technique to an Af-body setting, and it is demonstrated 
that this is indeed the natural way to do this. Alternatively, in a more general context, 
this formalism allows the study of a system which gradually loses particles to some 
environment due to the non-Hermitian part of // e ff . 

We comment that the absorption process does not affect the dynamics on its 
interior. However, removing a particle may still affect the other particles via their 
interaction. Hence, in order to obtain physically meaningful results, any dependence of 
the absorber must be eliminated by placing the absorber sufficiently far away from the 
interaction region - as is the case for any implementation featuring absorbing boundary 
conditions. 

In the following section, section [2l the formalism will be formulated for identical 
fermions exposed to a complex absorbing potential. In section |3] two examples 
featuring two identical fermions in one dimension are given. Finally, in section HI 
conclusions are drawn and a few future perspectives are outlined. 

2. Fock space description 

The natural setting for describing a system of a variable number of fermions is Fock 
space, the direct sum of all n-fermion spaces we wish to describe at most N 

particles, it suffices to consider 

/V 

^ = ©JC (8) 

n=0 

An arbitrary n-fermion state E can be written 

|¥) = I U^{xvx n )^{x x )---^{x n )\-), (9) 
n\ J 

where the field operator ij)\x) creates a particle at position x. The field operators obey 
the usual anti-commutator relations 

{^{x),^{x')} = 0, {T/>(x),t/>V)} = S(x-x'). (10) 

Here, x may denote all degrees of freedom associated with a particle, including 
eigenspin. Moreover, *P(xi ■■■x n ) is the anti-symmetric local wave function of the 
77-particle system, and | — ) G J4?o is the vacuum state, containing zero particles. 

A system of n identical fermions interacting with at most two-body forces has the 
Hamiltonian 

n n 

H = Y,h(x i )+ J £u(x„x j ). (11) 

j=l j<i 

Here, h{xi) (resp. u(xi,Xj)) is a one-body (two-body) operator acting on the degrees- 
of-freedom associated with particle i (and j). Using field operators, the Hamiltonian 
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can be written 

H= J dx^(x)h(x)^(x) + (12) 

- J dxdx' il>\x)if>\x')u(x,x')il>(x')if>(x), 

the point here being that H given on this form is independent of the number of particles 
in the system. The exact expression is irrelevant at this point. The CAP, which is 
diagonal in x, is conveniently introduced as 

- iT = -ij dx^(x)T(x)tf)(x). (13) 

By comparing ©, © and (fT3l) with the von Neumann equation ©, we see that 
the Lindblad operators, for a diagonal representation of the Lindbladian, must fulfil 

T = J dxA\x)A(x) (14) 

in order to reproduce the anti Hermitian "interaction" leading to absorption. Here 
we have allowed for an integral instead of a sum in ©. Furthermore, the last term 
of the Lindbladian, which is absent in the von Neumann equation, —2A n pA\ — > 
— 2A(x)pA^ (x), should map an Af-particle system into an (N — l)-particle system. 
Hence, A(x) should map J% into J4?n-u which means that the Lindblad operator 
A(x) must be on the form 

A(x) = J dxa(x,x')ip(x / ). (15) 

The simplest choice that satisfies (fl4)) is the diagonal form 



A(x) = y/T(x)il>(x). (16) 

We will return to the justification of why this is the proper way to define the Lindblad 
operators shortly. 

With (fT6l) . the Lindbladian CO), may immediately be written 

@(p) = iTp + ipT-2i J dxY{x)^{x)pil)\x), (17) 

and the master equation © becomes 

ih—p= [H,p] - i{T,p} + 2i j dxT(x)^(x)p^ (x). (18) 

This is our fundamental dynamical formulation of particle loss due to a CAP. 

Let us consider the master equation (fT8~l) in some detail. We may partition the 
density matrix p into blocks, viz, 

N N 

P = £ £ Pn,m, (19) 
7?=0 m=0 

where p W)OT = P n pP m , with P n being the orthogonal projector onto M* n . Each block 
evolves according to 

ih—p mt „ = [H,p mtn ] - i{T,p nhn } + (20) 



2iJ dxT(x)i}>(x)p m+ ^ n+ i^ '(x) 
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showing that the flow of p n>m depends on that of p, 7 + 1 , m + 1 , but not the other way around. 
This is illustrated in figure [Q Also, notice that obeys the original von Neumann 
equation © as there are no couplings to other blocks in this case as particles do not 
enter the Af-particle system. 

We now return to the question of whether the definition (fT6l) of the Lindblad 
operators, which obviously is the simplest one, is the only adequate one. Indeed, 
during the transition from an iV-particle system to an (N — l)-particle system, only 
the unabsorbed part of the original system should be reproduced within the (N — 1)- 
particle system, and this leads to the above choice. To see this, consider a somewhat 
idealized example consisting of two non-interacting particles. Since they do not 
interact, their wave function is given by a Slater determinant at all times, *F = 
[a(xi,t)/3(x2,t) —J3(xi,t)a(x2,t)]/V2. The state a does not overlap with the absorber 
at any time, i.e. 



and /3 corresponds to an unbound particle travelling outwards. We suppose that as 
t approaches infinity, the particle in the state yS is completely absorbed, and our 
numerical representation of the final system should converge towards the pure one- 
particle state a. 

Indeed, the evolution dictated by (1201) follows this pattern. The block of the 
density matrix corresponding to the two-particle system, p2,i, remains a pure state 
- albeit with decreasing norm as f3 is absorbed - and the evolution of the one-particle 
block pu simplifies due to (f2TT) to 



Hence, the one-particle part of the system is always proportional to a pure a-state, 
and, since the trace of the entire system is conserved, this simply integrates to 
Pi,i(? — > oo) = |a)(or|, as it should. 

On the other hand, with a non-diagonal form of the Lindblad operators, c.f. (fT5l) . 
we would have contributions to pij of the form \/3}(a\ and its Hermitian adjoint in 
addition to the pure state contribution. These contributions clearly cannot be allowed, 
as pi i should be independent of the state {3 of the outgoing particle. Hence, we find 
that, up to an arbitrary phase factor, the definition (fT6l) is the proper way to define A (x). 

For a non-interacting iV-fermion system where N a particles are bound and Np 
particles are unbound, and with an initial state given by the Slater determinant \m) — 
\{ai,...,a^ a ,/3i,...,/3^}}, it may be verified by inspection that the source term for the 
N a -th block, pN a ,N a , is always proportional to \{a\,...,a Na }){{a\,...,a Na }\ given that 
none of the a ; -states overlap with Y. The intermediate density matrix blocks PN-n,N-n, 
where < n < Np, will approach zero as t — > °°, but transiently describe (mixed) states 
where n particles have left. 

Of course, in the more interesting context of interacting particles, the structure of 
the diagonal blocks of the density operator, p n>n , is more complex than in the special 
case of non-interacting particles. 



T(x)a(x, t) = 0, for all x and t. 



(21) 




(22) 
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2.1. Consequences and interpretation 

Typically, our starting point is a pure Af-particle state, 1^(0) ) (^(O) | , in which case (1201 ) 
reduces to the ordinary Schrodinger equation ©, which was our original formulation. 
Moreover, it is easy to see that p(t) will remain block diagonal for all t in this case, 
i.e. p nm = if n ^ m. But (1201) shows that p n ,n{t) in general is a mixed state, unlike 
PN,N{t) = 1^(0) (^(0 1- This is due to the information loss when admitting ignorance 
of the whereabouts of the removed particle. 

We stress that by construction, Tr[p(f)] = 1 for all t. Probability flows 
monotonically from Jtffa into J^, and the particle absorbed from pm,n is not present 
in p n>n , but is erased, and p„ 5 „ is a proper description of an n-fermion system. In this 
way, we see that the above construction is a natural generalization of the original non- 
Hermitian dynamics, which describes the classical removal of a particle, into one that 
also describes the remaining system in a consistent way. It is striking to notice that the 
CAP, — zT(jc), is already given as one of the terms in the Lindbladian, so that the Fock 
space formulation follows almost immediately. 

It is worthwhile to mention that the traces of the blocks p„ jM along the diagonal of 
p have very simple interpretations as the probability P(n) of having n particles in the 
system upon a measurement, i.e. 

P(n)=Tr n \p(t)]=Tr\p n>n (t)]. (23) 

In particular, P(N) = (W^) < 1, and P(0) = 1 -Y% =1 P(n). Hence, within 
this formalism, distinguishing between single, double etc. ionization of atoms and 
molecules is straightforward. 

For any observable A the expectation value is given by (A) = Tr[Ap] = Tr[Ap„] . 
For example, the expected number of particles in the system is given by 

r N 
(jf) = dx Tx[^{x)^{x)p\ = £ nP{n). (24) 

It may also be useful to consider conditional expectation values: 

(A), = (25) 

i.e., the expectation value of A given that the system is found in an n-particle state. 

Obviously, as particles are "removed" by the absorber, information is lost. This 
information loss may be quantified by the von Neumann entropy, S = — Tr[plogp], or 
the closely related notion of purity, defined by lfl4l 

<r = Tr(p 2 )<l. (26) 

Unity minus this quantity, 1 — ^, is a measure of the amount of mixedness, and the 
purity g is 1 for pure states only. Similar to conditional expectation values, c.f. (T231) . 
one may define conditional purity and von Neumann entropy as the corresponding 
quantity of the re-normalized block, i.e. p,;,«/Tr[p n n ]. 

Of course, for a full quantum mechanical description in terms of the complete 
(unabsorbed) Af-particle wave function, no information is lost. Indeed, the absorption 
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of particles is a semi-classical concept. In the full A^-fermion quantum system no 
such separation is possible. On the other hand, as the Schrodinger equation for the 
N particles is separable in the non-interacting or, and /3, states discussed above, giving 
a Slater determinant of the time evolved one-particle states as the full solution, the 
Lindblad equation is seen to correctly construct the A^-particle Slater determinant 
resulting from the removal of the Np outgoing particles from this Slater determinant. 
Thus, the Lindblad equation exactly encapsulates the approximate separation of non- 
interacting quantum systems far apart. 

3. Example: Two identical spin ^-particles in one dimension 

In the following we consider two fermions in one dimension with the one-body 
Hamiltonian h given by 

h 2 d 2 

h = - — -^ + V(x,t), (27) 
2m cur 

where V(x,t) is some external, possibly time-dependent, potential. The fermions 
interact via a potential U (jq — x 2 ) and the extension of the system is effectively reduced 
by imposing a CAP. With two interacting identical fermions, we are confined to the 
Hilbert space 

3e = je 2 ®M[®M. (28) 

We discretize this system using a uniform grid in the interval [0, x max ) containing 
N points {xjj^Q , xj = jh, with h = x max /N. The field operators 1ft (x) can be replaced 
by a finite number of creation operators c -, creating a particle at grid point Xj. These 
operators, which obey the usual fermion anti-commutator rules 

{cj,c k }=0, {c j ,cl} = 6j >k , (29) 

map discrete anti- symmetrized ^-function bases for Jtf r into bases for M' n ±\ . 
The Hamiltonian takes the form 

H = E V? c ; + \ E u (*i - Xj)clc]cjCi, (30) 

where hij are the matrix elements of the one-body Hamiltonian, which depend on 
the chosen discretization of the second derivative. We choose a typical spectral 
approximation using the discrete Fourier transform, which also imposes periodic 
boundary conditions. The CAP is similarly discretized as 

-iT=-iJ^r(xj)c}c Jt (31) 

where T(x) is a non-negative function which increases as one approaches the boundary 
of the interval [0,x max ) and is zero in most of the interior. 

Concerning the density operator p, we write p = + p\ + Po, with p n = 

P n pP n . Initially the system is prepared in a two-particle pure state localized inside 
the absorption-free part of the grid. The master equation for p%{t) is equivalent to the 
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usual Schrodinger equation © for ^(f)), while the master equation for p\ (t) acquires 
a source term, i.e., 

ihjflit) = [H,pi(t)] ~ {r,pi(0} + (32) 

2i^r( Xj )cjmt)){^(t)\c]. 

j 

It is natural to view discrete the one-particle wave function as a vector </q (xf) of 
length equal to the number of discrete degrees of freedom. Similarly, it is natural to 
view two-particle wave functions as anti-symmetric matrices ^2 (•*,/' Moreover, a 
one-particle density matrix becomes a Hermitian matrix p\ (x, x') . Using these notions, 
the master equation (l32l) can be compactly written as 

ihj t Pi(t) = [H,pi(t)]-{r, Pl (t)} + 2ip s (33) 
with p s =2/zt/4l>2, 

where the last term contains matrix-matrix products. The extra factor 2 in p$ originates 
from the fact that the matrix 1^2 relates to a (redundant) "basis" of direct product states. 

If the initial two-particle state is an eigenstate of the total spin and its projection, 
the source term is always proportional to a single one-particle spin state. Hence, the 
spin does not introduce any complication in the notation in this case. For parallel 
spins, the one particle spin has the same direction as the two-particle spins, and for 
spin projection Ms = the one-particle density operator has its spinor component 
given 

by (I t) + (-l) S \l))/V2, where S is the total spin eigen value. 
The equation forp = Po(t)\— )(— | becomes 

^Po(t) = l%nxj)c jPl c]. (34) 

In principle, it is not necessary to include this equation in our calculations, as po(t) 
can be calculated from the constraint Tr[p(f)] = 1. 



3.1. Collision in a Gaussian well 

We will now focus on an example in which we set h = m= 1 and place the electrons 
in a potential of Gaussian form, 

2a 1 

The particles interact via a regularized Coulomb interaction 

U{\x l -x 2 \) = i (36) 

'xi-x 2 ) 2 + 6j 



VW = -\/oexp(-^-? L )- (35) 



c 

For this problem we chose a CAP of power form: 

r« = c(£)", 07) 

£ = max{0,x T -x,x-(x meix -XT)}, 
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where n > 1 and xj is the distance from the edges at which the CAP is "turned on", 
see figure [2] 

The system may for instance serve as a model for a quantum dot which couples 
to the conduction band and has narrow confinement in two dimensions. 

Equation (1331) is integrated using a scheme of second order in the time step based 
on a standard split-step operator scheme lfT8l . It is instructive to consider a more 
general setting in order to introduce the time stepping scheme for the density matrix. 
Given a differential equation for an entity y(t) in a linear space on the form 

y(t)=L t y(t)+f(t), (38) 

where L t is a linear operator dependent on t and f(t) is a source term independent of 
y(t), we may integrate formally using standard time-ordering techniques to obtain 

y (t) = ^e^ ds y(0)+F(t) + 

■■L Sn F(s n )ds n ---dsu (39) 

with F(t) = jQf(s)ds. The source terms are not easily transformed using the time- 
ordering operator 3 '. Assuming that the case f(t) = can be integrated to p-th order in 
the time t using y(t) = %y(0), a p-th order method for the general case can be obtained 
by keeping p source terms and evaluating these with sufficiently high order quadrature. 
For example, using standard Strang splitting which gives a scheme of local error 0(t 3 ), 
or other schemes based on the Magnus expansion |fl9ll . we may approximate the term 
F(t) as F(t) ps t[f(0) +/(/)] /2 (using trapezoidal quadrature) and the n = 1 term as 
f L s F(s)ds ~ t 2 L t /2[f(0) + tf (0)/2\. Specifically, in our implementation we have 
used a second order scheme given by 

Pl (j) = e-W+ u -^ tl2 e- lTt e-^ +u -^ tl2 Pl (0) x 

e+ i(V + U +l T) t /2 e -iTt e +i(V+U+iT)t/2 + ^ (Q) + 

t 2 \ps(0)-i(Hp s (0)-p s (0)H 
where ps = 2h (^Tiffl + i/^ri/rt 




+ 0(t 3 ) (40) 



Here T is the kinetic energy. As this scheme is trace preserving to second order only, 
we have also solved (l34l) in order to check that our numerical time step is small enough 
to preserve the total trace. 

Figure [3] shows the evolution of the particle density for a system in which a 
fermion collides with another one. In this case the potential depth Vq = 4 and the 
width cr = 0.75, the interaction strength A = 5, and the "softening" 5c = 0.1. The 
CAP has the power n = 3, the strength C = 4 and xj = 5. The starting point is a 
spatially anti-symmetrized state (spin triplet) consisting of a particle trapped in the 
well and an incoming wave packet of Gaussian shape. The trapped part corresponds to 
a superposition of the ground and the first excited one-particle states in the well, and the 
incoming wave function has mean momentum &o = 2. It is seen that as absorption, both 
due to transmission and back-scattering, takes place, the two-particle density vanishes 
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and a one-particle density emerges. It is also seen that, apart from in the absorption 
region, the total particle density, obtained by adding the two and one-particle densities, 
compares well to the "true" particle density obtained from solving the Schrodinger 
equation without absorber on a larger grid. 

Due to the collision, there is a finite probability that both particles are absorbed. 
This is clearly seen in figure HI which shows how the total trace is distributed between 
the sub-spaces Jfy, Jt\ and J% as a function of time. In this particular case we have 
P(l;t — > °o) = 0.92 and P(0;t — > oo) = 0.077. Also shown are the expectation value 
of the particle number and the purity of the system. Purity is reduced in two ways. 
Firstly, it is reduced as the trace becomes distributed between the three sub-systems, 
and secondly because p\ is not a pure state within This is seen from the fact that 
conditional purity, q\, converges towards 0.6, i.e. less than unity (not shown explicitly 
in the figure). 



3.2. Laser ionization of a one -dimensional helium atom 

In this next example we expose our system to an electric pulse of type 

£(0=£ sin 2 (^)cosM- (41) 

With this time-dependent perturbation, the one-particle Hamiltonian (l27l) acquires a 
time-dependent term, which in the length gauge representation reads 

xE(t). (42) 

We have here set the charge of the particles (the electrons) to be —1. The static 
potential Vn felt by the particles is chosen to be a regularized Coulomb potential, 

V N (x) = , (43) 

whereas the interaction between them is still described by (l36l) . By choosing 8 2 N — 
j a.u. and 8c = 0.5735 a.u., the ground state energy and the first ionization threshold 
coincide with those of a true three-dimensional atom, i.e. the ground state energy is 
—2.904 a.u. and the ground state energy of "He + " is —2 a.u.. By "a.u." we mean 
atomic units, defined by choosing the Bohr radius, the electron mass, the elementary 
charge and fi as units for their respective quantities. The ground state of the system, 
which is a spin singlet state, is easily obtained by propagation in imaginary time. 

Figure |5] shows the evolution of the system exposed to a pulse of maximum field 
strength Eq = 5 a.u., central frequency w = 3.2 a.u. and a duration corresponding to 
three optical cycles. This central frequency corresponds to a photon energy which 
energetically allows for one photon double ionization. Rather than an absorber of 
power-form (1371) . we have here used a Manolopoulos-type absorber 0, which has the 
advantages that it is transmission free and dependent on one parameter only. Along 
with a figure showing how the partial traces, given by (1231) . evolve in time, and another 
one showing the time-dependence of the electric field, given by (|4T|) . we have included 
snapshots of the absolute values of the two-particle wave function ^2(^1^2) an d the 
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one-particle density matrix p\(x,x'). It is clearly seen that as absorption takes place in 
the two-particle sub-system, a one-particle density matrix emerges. Note that the lobes 
following the axes of the {x,x!) -coordinate system do not correspond to absorption of 
the second electron but rather loss of coherence within the one-particle sub-system. 
However, from the lower panel of figure [51 which clearly demonstrates how single 
ionization may be distinguished from double ionization, we see that there is a finite 
probability of absorbing both particles. In this case, specifically, the probability 
of ionizing only one electron is P(l;t — > °°) = 0.31, and the probability of double 
ionization is P(0; £ °°) = 0.034. 

4. Conclusion 

In conclusion, we have demonstrated how the Lindblad equation in Fock space can be 
used to generalize the notion of absorbing boundary conditions for Af-particle systems. 
Specifically, the remainder of the system is preserved as a particle is absorbed. With 
this formalism it may be possible to describe the dynamics of unbound systems which 
otherwise would require an unrealistically large numerical grid. 

As a consequence of this being a Markovian process, some coherence is lost, 
and the state after absorption is in general given by a mixed state rather than a wave 
function. Within Lindblad theory, the identification between the Lindblad operators 
and the creation and annihilation operators comes quite natural for complex absorbing 
potentials. We believe that the method outlined in this paper may be generalized to 
other kinds of non-Hermitian "interactions" as well. 

Since it is considerably more involved to solve a master equation rather than a 
Schrodinger equation, c.f. EOl . future work will aim to find lower rank methods for 
solving (1201) . Also, more sophisticated spatial approximations like sparse grids |[T9l 
may be utilized to simplify the treatment of more particles. 
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Figure 3. Colour online: The evolution of the particle density. The upper panel is the 
particle density within the two-particle system, the second one from above is obtained 
from the one-particle density operator, and the third panel is the total particle density. 
Also shown, lower panel, is the particle density in the same region obtained from a 
solution of the full two-particle problem on a larger grid without absorber. The initial 
state is an anti-symmetrized state consisting of a Gaussian wave travelling towards 
the right with mean momentum fco = 2 and a state corresponding to a particle initially 



Absorbing boundary conditions for dynamical many-body quantum systems 15 





---Two particles 
— One particle 

Vacuum 
— Total trace 

<n> 


\ / 
\ J 

I \ 
/ \ 
/ \ 



20 40 60 

Time 




20 40 60 



Time 



Figure 4. Colour online: The trace and the expectation value of the number operator 
(upper panel) and the purity of the system (lower panel). In both cases the separate 
contributions from the two, one and zero part of the total density operator have been 
shown. 



B, 42:125205,2009. 

[11] Kenneth C. Kulander. Multiphoton ionization of hydrogen: A time-dependent theory. Phys. Rev. 
A, 35(1):445^M7, 1987. 

[12] M Protopapas, CH Keitel, and PL Knight. Relativistic mass shift effects in adiabatic intense laser 

field stabilization of atoms. /. Phys. B, 29(16):L591-L598, 1996. 
[13] G. Lindblad. On the generators of quantum dynamical semigroups. Comm. Math. Phys., 48:119, 

1976. 

[14] M. Schlosshauer. Decoherence and the quantum-to-classical transition. Springer- Verlag, 2007. 
[15] A Sandulescu and H Scutaru. Open quantum-systems and the damping of collective modes in 

deep inelastic-collisions. Ann. Phys., 173(2):277-317, 1987. 
[16] Pawel Caban, Jakub Rembieliriski, Kordian A. Smolihski, and Zbigniew Walczak. Unstable 

particles as open quantum systems. Phys. Rev. A, 72(3):032106, 2005. 
[17] Reinhod A. Bertlmann, Walter Grimus, and Beatrix C. Hiesmayr. Unstable particles as open 

quantum systems. Phys. Rev. A, 73(5):054101, 2006. 
[18] M.D Feit, J. A Fleck Jr., and A Steiger. Solution of the schrodinger equation by a spectral method. 

J.Comput. Phys., 47(3):412 - 433, 1982. 
[19] C. Lubich. From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical 

Analysis. European Mathematical Society, 2008. 
[20] Jean Dalibard, Yvan Castin, and Klaus M0lmer. Wave-function approach to dissipative processes 

in quantum optics. Phys. Rev. Lett., 68(5):580-583, 1992. 



Absorbing boundary conditions for dynamical many -body quantum systems 



16 




20 30 
Time (a.u.) 



40 









• 


n 
* 


t 

— ♦ 




• 

• 








• 




Figure 5. Colour online: The panels show the evolution of a model one-dimensional 
helium atom exposed to a five cycle sine-squared electromagnetic pulse of strength 
Eq = 5 a.u. and central frequency <±> = 3.2 a.u. (upper panel). The middle panels show 
the absolute value of the wave function of the two-particle part ^2(^1 > x 2) (left) an d the 
absolute value of the density matrix of the one-particle part p\ (x,x r ) (right) at various 
instances. The time of each "snapshot" is indicated by a diamond in the upper panel. 
The lower panel shows the probability of finding two (dashed curve), one (full curve) 
and zero particles (dash-dotted curve) in the system in the same manner as in the 
upper panel of figure|4] The dotted curve shows the probability of the system being in 

ttif ini 



